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Abstract 

The  status  of  computer  simulations  of  electric  double  layers  is  briefly  summarized  and  a  road 
map  with  bottle  necks  for  solving  the  important  problems  in  the  atomic  scale  simulation  of 
interfacial  electrochemical  processes  is  proposed.  As  an  example  of  recent  activity,  efforts 
to  simulate  screening  in  electric  double  layers  are  described.  Molecular  dynamics  simulations 
on  systems  about  4  nm  thick,  containing  up  to  1600  water  molecules  and  NaCl  at  IM  to 
3M  concentrations,  are  shown  to  exhibit  the  main  components  of  electric  double  layers  at 
charged  metal  surfaces.  Ion  and  water  density  profiles  across  the  system  show:  a  bulk 
electrolyte  zone,  a  diffuse  layer  that  screens  the  charge  on  the  electrode  and  a  layer  of  ori¬ 
ented  water  on  the  surface. 
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The  status  of  computer  simulations  of  electric  double  layers  is  briefly 
summarized  and  a  road  map  with  bottle  necks  for  solving  the  important 
problems  in  the  atomic  scale  simulation  of  interfacial  electrochemical 
processes  is  proposed.  As  an  example  of  recent  activity,  efforts  to 
simulate  screening  in  electric  double  layers  are  described.  Molecular 
dynamics  simulations  on  systems  about  4  nm  thick,  containing  up  to 
1600  water  molecules  and  NaCl  at  IM  to  3M  concentrations,  are 
shown  to  exhibit  the  main  components  of  electric  double  layers  at 
charged  metal  surfaces.  Ion  and  water  density  profiles  across  the 
system  show:  a  bulk  electrolyte  zone,  a  diffuse  layer  that  screens  the 
charge  on  the  electrode  and  a  layer  of  oriented  water  on  the  surface. 


This  paper  describes  the  application  of  molecular  dynamics  to  chemical  processes  at 
the  interface  between  a  charged  metal  electrode  and  aqueous  electrolyte.  The  long 
range  goal  is  a  scheme  for  the  dynamics  of  chemical  reactions  on  surfaces  important 
in  the  electrochemical  technology  of  power  sources,  electroplating,  and  corrosion 
control.  The  paper  begins  with  a  summary  of  our  view  of  the  current  state  of  com¬ 
puter  simulation  applied  to  interfacial  electrochemistry.  Next  we  outline  a  map  on 
which  we  plot  a  mdimentary  road  towards  the  goal.  The  status  and  road  map  are 
accompanied  with  a  commentary  on  critical  problems  and  potential  bottle  necks.  To 
illustrate  progress  in  the  field  we  describe  our  simulations  of  screening  of  charged 
electrodes  by  aqueous  electrolytes  including  previously  unpublished  work.  Screening 
by  double  layers  is  an  important  physical  phenomenon  because  double  layers  are 
some  of  the  basic  organizations  adopted  in  electrochemical  and  biological  systems  to 
shield  electric  fields  arising  from  layers  or  arrays  of  charge  in  contiguous  structures. 
It  is  therefore  important  to  understand  their  properties  using  models  that  can  be  solved 
without  further  approximations.  The  structure  of  the  aqueous  part  of  the  double  layer 
is  given  in  terms  of  time  independent  water  and  ion  probability  distribution 
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functions  averaged  parallel  to  the  metal  surface.  Electric  fields  and  potentials  are 
calculated  from  the  microscopic  charge  density  profile.  These  calculations  provide 
a  consistent  microscopic  picture  of  ions  and  water  throughout  a  double  layer  includ¬ 
ing  the  species  next  to  the  charged  surface  (inner  layer),  in  the  'diffuse  layer*  (also 
called  the  screening  layer)  and  in  the  bulk  zone.  The  effect  of  finite  sized  ions  and 
water  are  clearly  evident,  as  is  the  effect  of  the  electric  field  on  the  orientation  of 
surface  water  molecules. 

Status  of  Molecular  Dynamics  Simulations 

Figure  1  (top)  is  a  sketch  of  the  electric  potential  using  Gouy-Qiapman-Stem 
theory[l]  in  which  the  diffuse  layer  ions  are  treated  as  point  charges,  the  water  as  a 
a  dielectric  continuum,  and  OHP  and  IHP  (outer  and  inner  Helmholtz  planes)  are 
introduced  to  marie  the  distance  of  closest  approach  of  strongly  hydrated  ions  and 
contact  adsorbing  ions  respectively,  molecules  and  ions  for  times  as  long  as  several 
nanoseconds.  Research  is  not  limited  by  computer  power  but  by  the  availability  of 
correct  theories  Figure  1  (bottom)  shows  a  molecular  scale  cartoon  of  ions  and  water 
near  a  flat  charged  metal  surface.  Ideas  embodied  by  pictures  like  this  together  with 
Gouy-Chapman-Stem  theory  and  the  thermodynamic  theory  of  surface  excess  quan¬ 
tities  have  been  used  to  analyze  and  interpret  experimental  electrochemical 
data[2-5].  The  advent  of  rise  based  work  stations  allows  the  testing  of  atomic  scale 
models  of  interfaces  with  hundreds  and  even  thousands  of  molecules.  Monte-Carlo 
and  molecular  dynamics  computer  simulations  of  ions  and  water  molecules 


Figure  1.  Schematic  diagram  (top)  of  electric  potential  across  the  double  layer  and 
(bottom)  cartoon  of  the  structure  of  water  and  ions  next  to  a  flat  charged  electrode. 


interacting  through  various  named  potentials  are  routinely  performed  for  a  few  thou¬ 
sand  water  of  interfacial  chemical  processes  and  good  algorithms  for  the  efficient 
calculation  of  sums  of  long  range  electrostatic  interactions.  We  comment  first  on 
electrostatics  and  boundary  conditions,  and  then  move  on  to  other  scientific  issues. 

Electrostatics.  Three  dimensional  sums  of  electrostatic  interactions  (eg.,  ion-ion, 
dipole-dipole)  are  conditionally  convergent.  Available  algorithms  are  getting  faster 
and  more  useful  all  the  time.  However  in  studying  the  physics  of  electrified  inter¬ 
faces  it  is  essential  that  the  long  range  part  of  the  electrostatic  interaction  be  computed 
without  truncation  in  a  manner  consistent  with  the  boundary  conditions.  For  example 
on  a  metal  electrode  Maxwell's  equations  connect  electric  field,  potential  with  the 
total  charge  distribution  (metal  surface  charge  and  charge  on  molecules  and  ions). 
For  some  geometries  (eg.,  planar)  the  electric  field  of  surface  charge  can  be  calculated 
by  the  method  of  images.  In  3D  systems  if  the  electrostatic  interaction  is  truncated 
the  long  range  correlations  responsible  for  dielectric  polarization  will  be  wrongly 
computed  possibly  causing  like  charged  systems  to  attract.  Parallel  metal  siufaces 
have  an  infinite  set  of  multiple  electrostatic  images  that  have  to  be  summed  in 
plane-wise  fashion  to  get  the  physics  right[6,  7].  The  crystal  optics  based  methods 
of  Ewald  and  Komfeld  are  the  simplest  for  calculating  electrostatic  fields[8].  The 
algorithm  woiks  for  all  space  group  symmetries.  At  best  diis  method  is  order  iV  *  in 
the  number  of  charges  N.  There  are  2D  summation  methods  that  are  faster[7,  9]. 
In  our  simulations  we  use  the  order  N  fast  multipole  method  (finm)  developed  by 
Greengard  and  Rokhlin  [10-13].  It  is  a  very  useful  method  for  large  electrochemical 
simulations  where  a  variety  of  boundary  conditions  (periodic,  Dirichlet,  Neumann  or 
mixed  boundaries)  are  encountered.  It  can  also  be  adapted  so  that  regions  of  low 
charge  density  are  not  subdivided  when  the  charge  count  falls  below  a  specified  in- 
teger[12].  It  is  restricted  to  cubic  shaped  cells.  The  fmm  is  faster  than  Ewald  for 
systems  exceeding  a  few  hundred  charges[14].  Particle-mesh  methods  have  been 
extensively  used  for  long  range  r*'  potential  problems[15].  For  most  systems  P^M 
is  faster  than  fmm  and  can  be  used  with  orthorhombic  simulation  cells[16].  Recently 
we  showed  how  to  correctly  calculate  the  fields  from  charge  distributions.  The  effect 
of  averaging  over  space  regions  spanning  the  size  of  a  water  molecule  was 
explored[17,  18].  Though  most  experiments  are  done  at  constant  potential,  there  are 
very  few  simulations  at  constant  potential[19,  20]. 

Molecule-molecule  and  molecule-metal  potentials.  There  are  continuing  improve¬ 
ments  in  molecule-molecule  potentials.  High  quality  efforts  are  directed  at  improving 
the  interaction,  including  polarizabilies  and  getting  away  from  vanilla  Lennard-Jones 
forms[21,  22].  There  are  also  potentials  that  include  three  body  terms 
explicitly[23].  Possibly  the  best  atom-metal  potential  is  due  to  Barker[24].  The 
Barker  potential  for  Xe/Pt(lll)  is  an  excellent  fit  to  a  large  body  of  experimental 
data.  There  have  been  numerous  quantum  chemistry  studies  of  ions  and  water  on 
metal  clusters  some  with  applied  fields  and  others  on  charged  clusters  to  imitate  the 
electrochemical  environment  (for  water  references  see  Zhu[25],  for  ions  on  clusters 
see  Pacchioni[26]).  Unfortunately  only  a  few  have  been  parameterized  into  pair 
forms  useable  in  an  MD  code[25,  27-29].  For  electrochemically  interesting  metals 
like  Pt  there  are  many  compute  intensive  problems  associated  with  cluster  based  po- 


tentials  of  third  row  transition  metals.  Not  the  least  are  the  relativistic  orbital  con¬ 
tractions  and  spin  orbit  effects.  It  would  be  better  to  focus  on  the  sp  metals  like  Cu 
or  Ag.  Even  Au  would  be  easier  to  work  with  than  Pt  since  its  d  shell  is  more  tightly 
bound.  Recently  several  publications  have  reported  cluster  calculations  for  water  and 
ions  on  Hg[30,  31]  including  a  parameterization  to  give  a  pair  potential  for  MD 
studies. 

Dynamics.  There  are  calculations  in  which  the  metal  is  modeled  as  an  Einstein  solid 
with  harmonic  vibrations[32].  When  surface  molecules  and  ions  are  strongly 
adsorbed  molecular  dynamics  becomes  an  inefficient  way  to  study  surface  processes 
due  to  the  slow  exchange  between  surface  and  solution.  In  this  case  it  is  possible  to 
use  umbrella  sampling  to  compute  distribution  profiles[33,  34].  Recently  the  idea 
underlying  Car-Parrinello  was  used  for  macroion  dynamics[35,  36]  in  which  the 
solvent  surrounding  charged  macroions  is  treated  as  a  continuum  in  a  self  consistent 
scheme  for  the  potential  controlling  ion  dynamics.  Dynamical  corrections  from  the 
solvent  can  be  added.  There  is  a  need  to  develop  statistical  methods  to  treat  the  dy¬ 
namics  of  complex  objects  that  evolve  on  several  different  time  scales. 

Interfacial  electron  transfer.  There  have  been  several  studies  of  electron  transfer 
reactions[19,  37]and  the  connection  with  Marcus's  theoiy[38].  It  may  be  possible  to 
use  a  Car-Parrinello  like  scheme  on  that  part  of  the  system  directly  affecting  the 
electron  transfer.  There  has  also  been  very  interesting  studies  of  the  ferro-ferri  redox 
couple  in  solution[39,  40]  that  address  many  issues  related  to  electron  transfer  from 
an  electrode  to  a  hydrated  ion.  Slow  processes  can  be  treated  by  transition  state 
methods  like  the  ones  used  in  solid  state  ionic  conductivity [41]. 

Ensembles.  One  goal  of  simulations  is  the  calculation  of  experimental  quantities. 
The  most  challenging  are  the  Gibbs  free  energies  of  adsorption.  Currently  there  is 
no  proven  scheme  for  constant  chemical  potential  simulations  of  electrolyte 
adsorption  on  metals.  It  seems  possible  to  develop  Andersen's  method[42]  for 
(N,P,T)  ensembles.  Another  possibility  is  an  analog  of  the  Gibbs  ensemble[43]  for 
electrolytes  between  plates  with  a  bulk  sample  forming  the  second  phase.  In  an  in¬ 
teresting  recent  development  the  grand  canonical  Monte  Carlo  method  was  used  for 
atomic  fluid  mixtures  in  a  slit  pore[44,  45]. 

Toward  a  Road  Map 

The  intention  is  modest  though  the  title  of  this  section  sounds  pretentious.  Given  the 
current  state  of  theory  and  simulation  can  we  identity  a  path  that  will  eventually  lead 
to  practical  computation  of  reactions  important  to  electrochemical  technologies.  To 
simplify  the  discussion  we  focus  on  technologies  connected  with:  electric  power 
sources,  electroplating,  and  corrosion  control.  The  key  science  problems  are  the 
following:  deposition  and  dissolution  of  metal,  formation  of  oxide  layers,  electron 
transfer  from  electrodes  to  ions,  and  charge  migration  in  complex  fluid  phases.  This 
is  a  broad  range  of  processes.  Let  us  first  examine  two  existing  methods:  ab  initio 
molecular  dynamics  and  dynamics  using  potential  energy  surfaces. 


Ab  initio  molecular  dynamics.  Chemical  reactions  involve  the  reorganization  of 
electrons  about  the  nucleii  involved  in  the  bond  changes.  The  ab  initio  molecular 
dynamics  scheme  developed  by  Car  and  Parrinello[35]  permits  an  accurate  de¬ 
scription  of  both  electronic  and  nuclear  rearrangements  that  occur  during  a  reaction. 
The  penalty  for  including  electronic  coordinates  explicitly  in  an  electrochemical 
simulation  is  the  restriction  to  relatively  few  water  molecules.  Liquid  water  and 
proton  transfer  have  been  studied[46,  47].  The  computational  problem  is  immense 
so  that  at  present  the  study  of  himdreds  of  water  molecules  takes  too  long.  This 
number  is  quite  enough  for  studying  H-bonding  and  dissociation  and  the  dynamics 
of  the  hydration  of  ions  but  is  insufficient  to  deal  with  double  layer  stmcture  or  re¬ 
actions  of  hydrated  ions  with  charged  metals. 

Potential  energy  surfaces.  For  systems  where  some  or  all  of  the  dynamics  can  be 
described  by  a  potential  energy  surface  (PES)  it  is  possible  to  avoid  solving  the 
electronic  Schrodinger  equation,  and  use  instead  a  PES  parameterized  with  exper¬ 
imental  data.  Several  cases  already  exist.  First  there  is  the  well  known  5/^  reaction 

cr  +  //sC-a  ->  C/-C//3  +  cr.  [i] 

The  molecular  dynamics  of  this  reaction  has  been  studied  in  the  gas  and  solution 
phases.  Another  example  is  the  Brenner  bond  order  potential[48].  used  to  describe 
the  dynamics  of  homolytic  bond  fission  and  formation  in  carbon  and  hydrogen  con¬ 
taining  systems.  These  examples  involve  the  making  or  breaking  of  two  electron 
bonds  between  low  Z  atoms.  Explicit  degrees  of  freedom  for  the  electrons  are 
avoided  though  the  use  of  a  PES.  The  electrochemical  problems  are  harder.  They 
involve  charged  metal  surfaces,  charge  transfer  in  polar  environment  and  dynamics 
on  several  time  scales.  In  surface  UHV  adsorbate  studies  considerable  progress  has 
been  made  understanding  chemical  dissociation  reactions  and  physical  sputtering  us¬ 
ing  a  combination  of  LEPS  potential  for  the  diatomic  on  the  metal  with  a  many  body 
embedded  atom  potential  for  the  metal. 

Road  map.  Consider  first  the  deposition  and  dissolution  of  metal  ions  on  charged 
surfaces.  If  systems  can  be  identified  where  adiabatic  potentials  could  be  used  to 
describe  the  nuclear  motion  then  by  analogy  with  the  experience  gained  with  the 
dissociation  of  diatomics  on  metal  surfaces[49]  and  with  carbon-hydrogen 
chemistry[50]  it  should  be  possible  to  describe  how  metal  ions  adsorb  on  metals  and 
conversely  how  metals  dissolve.  Once  this  is  achieved  it  should  be  possible  to  for¬ 
mulate  schemes  to  describe  the  formation  of  insoluble  metal  oxide  layers.  This  would 
also  require  a  model  for  water  decomposition  on  metals.  The  central  force  model 
which  permits  water  dissociation[51,  52]  would  have  to  be  developed  for  the  surface 
environment.  The  third  process  involves  understanding  how  electrons  cross  from  a 
metal  to  ions  or  organic  molecules  and  initiate  chemical  reactions.  It  is  possible  that 
some  of  these  reactions  could  be  studied  by  treating  the  electron  quantum  mechan¬ 
ically  moving  in  a  potential  defined  by  classical  mechanics  motion  of  molecules  and 
ions  in  the  double  layer.  Solvated  electrons  in  liquids  are  studied  this  way.  If  not 
then  hybrid  Car-Parrinello  schemes  would  have  to  be  developed.  The  fourth  class 
we  have  picked  is  important  for  ion  mobilities  in  lithium  ion  batteries.  It  may  be 


possible  to  study  these  systems  using  ideas  borrowed  from  small  polaron  motion  in 
solids,  using  transition  state  dynamic  methods  of  Bennett[4l],  or  using  recent  theories 
of  macroion  motions[36].  Though  we  are  presently  a  long  way  from  providing 
technologically  useful  information,  we  believe  that  being  able  to  model  aspects  of 
these  processes  would  provide  atomic  scale  insight  as  to  how  these  technologies 
work. 

Details  of  the  Model 

Screening  is  treated  using  the  immersed  electrode  model[14,  53,  54].  A  layer  of 
electrolyte  with  an  excess  of  positive  ions  is  confined  by  a  semi-infinite  metal  on  the 
rhs  and  a  restraining  non-metallic  wall  on  the  Ihs.  There  is  no  external  electric  field. 
The  metal  is  grounded  to  zero  electric  potential.  The  charge  on  the  electrode  equals 
the  image  charge  of  the  excess  positive  ions  (-4e  or  -8e)  so  that  the  metal  charge  is 
coiiq)letely  screened  by  the  ions  in  solution.  The  difference  between  these  and  earlier 
simulations  [14,  53-56]  is  size,  here  there  are  enough  ions  and  water  to  form  a  bulk 
electrolyte  region  where  the  solution  is  locally  neutral.  Because  there  is  only  one 
metal  surface  (shown  on  the  rhs  in  all  the  figures)  there  are  no  multiple  electrostatic 
images  in  the  direction  perpendicular  to  the  metal[6,  7].  The  electric  field  and  po¬ 
tentials  are  calculated  using  previously  developed  methods[17,  18]. 

This  work  uses  the  SPCE  water  model[57]  (three  charged  mass  points,  qjj  = 
0.4238e,  bond  angle  109.5®,  OH  bond  length  0.1  nm,  Lennard-Jones  sphere  with  ra¬ 
dius  cr  =  0.317  mn  and  well  depth  s=  0.650  kJ/mole)  and  associated  parameters  for 
NaCl[58].  The  coordinate  origin  is  at  the  center  of  the  cell,  and  the  axes  are  per¬ 
pendicular  (z)  and  parallel  (x,y)  to  the  metal  surface.  The  simulation  cell  has  edge 
length  L=3.724  nm,  the  flat  metal  plane  is  at  z  =  1.862  nm  and  the  flat  restraining 
wall  at  z  =  -1.862  nm.  The  cell  contains  upto  1600  water  molecules  and  the  ion 
concentrations  are  approximately  IM,  2M  or  3M  NaCl.  The  metal  charge  density  is 
either  -0.046  Cm'^  (-4eL"^)  or  -0.092  Cm'^  (-8eL‘^).  Contact  adsorption  of  ions  is 
minimized.  The  cation  has  a  smaller  radius  than  the  anion,  so  its  hydration  shell  is 
strongly  bound  making  it  difficult  for  it  to  contact  adsorb.  The  negative  metal  charge 
makes  it  energetically  unfavorable  for  the  anions  Cl"  to  contact  adsorb. 

Two  potentials  are  used  to  describe  the  interaction  of  water  and  ions  with  the 
metal.  A  9-3  potential  is  used  to  for  the  Pauli  repulsion  and  the  attractive  dispersive 
interactions  between  molecules  or  ions  and  the  metal.  The  interaction  between  a 
charge  on  an  ion  or  water  with  the  conduction  electrons  of  the  metal  is  modeled  with 
a  classical  electrostatic  image  potential.  The  position  of  image  plane  and  origin  plane 
(same  as  the  plane  through  nuclei  of  the  surface)  of  the  9-3  potential  was  taken  to 
be  coincident.  In  real  materials  the  image  and  nuclear  planes  are  not  coincident. 
This  is  not  important  in  the  simulations  because  the  thickness  of  the  repulsive  wall 
is  large  (ca.  0.247  nm).  The  9-3  atom-surface  wall  parameters  describing  interaction 
with  nonconduction  electrons  were  chosen  to  be  the  same  as  that  used  by  Lee 
etal[59],  A=17.447xl0'^  kJ(nm)Vmole  and  B=76. 144x10’^  kJ(nm)Vmole  for  the  O 
atom.  The  A  and  B  parameters  for  H  were  set  to  zero.  The  potential  corresponding 
to  these  parameters  describe  a  graphite-like  surface.  A  useful  reference  point  in  the 
wall  potential  is  at  z^=1.615  nm  where  the  9-3  wall  potential  changes  sign.  Each 
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simulation  was  run  for  about  a  nanosecond,  and  the  instantaneous  positions  of  all  the 
atoms  recorded  every  picosecond.  The  first  100  ps  were  discarded  as  anneal  time. 
The  density  probability  functions  p(z)  were  constmcted  by  binning  the  configurations 
in  bins  (width  L/800)  parallel  to  the  metal  surface  to  give  functions  of  z  only. 

Screening  of  Charged  Metal  Electrodes  in  SPCE  Electrolyte 

We  begin  with  remarks  on  the  screening  of  charged  surfaces  by  aqueous  electrolytes. 
At  high  salt  concentrations  the  region  with  excess  ionic  charge  is  microscopically 
small.  A  rough  estimate  of  the  thickness  of  the  zone  in  which  screening  occurs,  v^id 
for  dilute  solutions  (<  O.IM),  is  given  by  the  inverse  of  the  Debye-Huckel  screening 
constant[60]  _ _ 

d=^j(ekTI(SKe\v^))  .  [2] 

Here  s  is  the  macroscopic  dielectric  constant  of  the  solvent  (ca.  80  for  water),  v  the 
valence  of  the  ion  (one  in  this  paper),  and  is  the  bulk  concentration  of  the  ions. 
Typical  values  of  d  are:  3  nm  in  O.OIM,  1  nm  in  O.IM,  0.3  nm  in  l.OM,  and  0.2  nm 
in  3M  NaCl  solutions.  Obviously  in  high  salt  concentrations  die  screening  should 
be  more  efficient  and  the  screening  length  smaller,  but  since  the  number  for  IM  NaCl 
is  the  diameter  of  a  water  molecule,  and  the  number  for  3M  NaCl  is  even  smaller, 
these  numbers  mean  nothing  widiout  some  additional  atomic  scale  in  sight  like  that 
provided  by  MD  simulations.  At  high  concentrations  there  are  many  problems  with 
simple  Gouy-Chapman  theoty[61-63]  and  many  modifications  have  been 
proposed[60].  There  are  three  main  problems:  the  dielectric  constant  of  water  in  a 
high  surface  field,  the  lower  length  scale  due  to  the  finite  molecular  size,  and  corre¬ 
lated  motion  amongst  ions  and  water.  For  example  there  is  no  reason  to  believe  that 
the  value  of  s  is  80  near  a  surface  or  in  a  field  high  enough  for  dielectric  saturation 
to  occur.  This  aspect  has  been  discussed  many  times  in  the  electrochemical 
literature[64].  Though  we  take  the  Debye  length  for  concentrated  solutions  as  an 
rough  measure  of  double  layer  thickness,  we  should  be  very  cautious  when  d  ap¬ 
proaches  the  size  of  a  water  molecule. 

Figure  2  shows  the  distribution  of  the  ions  for  three  separate  calculations  with 
concentrations  IM,  2M  and  3M  NaCl.  The  charge  on  the  electrode  was  -4e  for  IM 
and  2M  and  -8e  for  3M.  The  temperature  was  30°C  for  IM  and  2M,  and  100°C  for 
the  3M  NaCl  solution.  Note  that  there  is  no  significant  contact  adsorption.  All  the 
peaks  in  the  ion  distributions  occur  away  from  the  position  of  closest  approach 
z  =1.612  nm  to  the  metal.  This  is  precisely  the  situation  we  contrived  for  the  reasons 
gTven  before.  The  ion  concentrations  are  approximately  the  same  for  |z|  <  1.0  nm. 
This  identifies  the  region  of  the  system  with  bulk  electrolyte  properties.  For  \z\  <  0.5 
nm  the  two  ion  profiles  are  the  same  within  10  to  20%,  for  all  the  simulations.  The 
bulk  region  is  smaller  for  IM  NaCl  because  the  screening  layer  is  thicker.  The  in¬ 
tegrated  ion  densities  are  monotonically  increasing  curves  in  Figure  2.  They  provide 
a  rough  measure  of  overall  charge  neutrality  from  the  restraining  wall  to  the  metal. 
The  vertical  arrows  indicate  roughly  where  the  ion  charge  densities  start  to  diverge 
because  of  screening.  Also  shown  in  Figure  2  are  the  results  of  a  calculation  of  ion 
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Figure  2.  Ion  probability  density  distribution  profiles  for  NaCl  solution.  IM  and 
2M  solution  at  30°C  and  -4e  electrode  charge,  3M  solution  at  100°C  and  -8e  electrode 
charge.  Also  shown;  integrated  ion  densities  and  Gouy-Chapman  ion  densities. 

densities  using  simple  Gouy-Chapman  theory.  For  optimal  comparison  the  electrode 
was  assume  to  start  at  z=1.612  nm.  These  superimposed  curves  show  how  lack  of 
atomic  scale  structure  limits  the  application  of  Gouy-Chapman  theory.  This  is 
graphic  evidence  that  there  are  important  details  in  MD  simulations  due  to  hydration, 
molecular  size,  and  water  layering  near  charged  surfaces. 

In  each  case  the  width  of  the  screening  layer  is  too  small  to  justify  the  de¬ 
scription  as  'diffuse'  (though  we  will  continue  to  use  this  term).  To  estimate  the 
diffuse  layer  width  we  measure  from  the  z=1.615  nm  (where  wall  potential  changes 
sign)  to  the  point  where  the  difference  in  integrated  ion  densities  is  e  '  of  the  metal 
charge  (8e  for  3M,  4e  for  IM  and  2M).  For  the  solutions  we  find:  0.5  nm  (0.31  nm) 


for  IM,  0.4  nm  (0.24  nm)  for  2M,  0.2  nm  (0.18  nm)  for  3M.  The  Debye  shielding 
lengths’ calculated  using  Eqn(2)  are  in  parentheses.  Made  with  these  favorable  as¬ 
sumptions  the  agreement  is  remarkably  close.  It  could  be  changed  immediately  if 
we  used  e  =  6,  a  value  more  appropriate  for  a  zone  in  which  the  dielectric  properties 
of  water  are  at  saturation  values[64]. 

Looking  at  the  fine  structure  in  the  density  profiles  we  see  that  on  the  metal  side 
all  the  chloride  distributions  have  weak  features  at  ca.  1.2  nm  and  1.4  nm.  Both 
appear  to  be  associated  with  peaks  in  the  cation  distribution  and  may  therefore  be 
due  to  contact  pairs  or  solvent  separated  pairs.  Of  course  the  SPCE  model  for  water 
was  not  designed  with  high  salt  concentrations  in  mind,  so  the  ion  pairs  may  be  more 
a  feature  of  the  model  and  not  nature.  Correlation  between  ions  at  high  salt  con¬ 
centrations  alter  the  distribution  near  the  charged  surface.  Note  that  there  is  no  as¬ 
sociation  of  oppositely  charged  ion  peaks  at  the  left  wall  in  Figure  2. 

IM  NaCl  Solution 

Figures  3,  4  and  5  show  the  results  of  an  MD  simulation  using  IM  NaCl  solution. 
In  this  simulation  the  cell  contains  32  Na^  ions,  28  Cl"  ions,  and  1576  water  mole¬ 
cules  at  30°C.  The  electrostatic  charge  on  the  electrode  surface  due  to  the  difference 
in  number  of  positive  over  negative  ions  is  -41e|  or  -0.046  Cm'^.  The  top  panel  in 
Figure  3  shows  the  probability  density  profiles  for  the  water  proton,  water  mass 
center,  Na"^  ions  and  Cl’  ions.  Both  ion  distributions  have  been  smoothed  to  permit 
clearer  identification  of  variations  in  density  with  position.  We  saw  in  Figure  2  the 
near  coincidence  of  the  integrals  of  the  ion  density  for  z  <  0.7  nm  which  shows  that 
the  electrolyte  is  approximately  charge  neutral  before  this  point.  For  z  >  0.7  nm  the 
integrated  densities  systematically  diverge  as  expected  for  a  transition  from  the  lo¬ 
cally  neutral  ’bulk’  electrolyte  into  the  'diffuse'  part  of  the  electric  double  layer.  The 
Na"^  ion  distribution  shows  well  defined  stmeture  in  the  form  of  a  broad  peak  at  ca. 
1.1  nm,  and  a  sharp  peak  at  1.4  nm.  The  water  and  proton  distributions  appear  flat 
for  \z\  <  0.8  nm.  On  the  metal  side  the  water  probability  distribution  has  peaks  at 
0.9  nm,  1.2  nm,  and  a  strong  asymmetric  peak  at  1.6  nm.  This  latter  feature  appears 
to  be  composite  being  the  superposition  of  a  broad  feature  at  1.5  nm  and  a  narrow 
peak  at  1.6  nm.  The  peaks  in  the  proton  distribution  are  at  0.9  nm,  1.2  nm,  1.55  nm, 
and  1.7  nm.  The  last  peak  at  ca.  1.7  nm  comes  from  protons  in  water  OH  bonds 
pointing  at  the  metal. 

The  bottom  panel  in  Figure  3  shows  the  potential  calculated  using  the  atom 
method,  and  the  components  of  the  potential  calculated  by  the  molecule 
method[17,  18].  The  contact  potential  is  about  -2V,  and  the  potential  in  the  'bulk¬ 
like'  zone  comes  from  the  water  quadrupole.  The  monopole  curve  is  from  the  ionic 
charge.  If  the  system  were  truly  neutral  then  the  monopole  curve  would  be  flat  and 
zero  all  the  way  to  the  beginning  of  the  diffuse  layer.  The  transition  to  monotonic 
decrease  starting  near  z  =  0.7  nm  is  another  indicator  of  where  the  diffuse  layer  be¬ 
gins  in  this  simulation.  Note  that  combined  monopole  and  dipole  potential  curve 
m+d  shows  that  the  dipole  potential  almost  completely  compensates  the  monopole. 
Adding  the  quadrupole  to  m+d  shifts  the  core  region  downwards  by  3V  and  brings 
the  molecular  calculation  of  potential  into  fair  correspondence  with  the  atom  method 


Potential  V  /  5.22  Density  profile  p  /  nm  ^ 


-1.6  -.8  0  .8  1.6 

Distance  z  /  nm 


Figure  3.  Top.  Probability  density  distribution  profiles  for  IM  NaCl  solution  at 
30°C  and  -4e  electrode  charge.  Bottom.  Total  electric  potential  and  component 
multipole  potentials.  Total  potential  by  atom  method. 
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re  4.  Detail  of  individual  electric  potentials  for  z  >  0.8  nm..  IM  NaCl  solution 
i°C  and  -4e  electrode  charge. 
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Figure  5.  Microscopic  electric  charge  density  and  electric  field  IM  NaCl  solution. 


of  calculation.  Adding  the  octopole  improves  the  agreement  at  the  walls.  The  mo¬ 
lecular  method  has  larger  extrema  near  the  surface  compared  to  the  atom  method. 
The  reason  for  this  is  the  need  to  include  many  high  order  multipoles  in  the  molecule 
method.  However  the  contact  potential  is  the  same  in  each  case  since  it  depends  only 
on  m  and  d[65]. 

Figure  4  shows  a  detail  of  the  the  potentials  shown  in  Figure  3  (bottom)  for  the 
region  z  >  0.8  nm.  We  note  that  the  monopole  and  dipole  determine  the  potential 
at  contact.  The  quadrupole  potential  is  very  important  in  bulk  and  at  the  surface. 
The  octopole  potential  is  important  near  the  surface.  In  the  future  as  water  models 
improve  it  will  be  very  important  to  include  these  terms.  Obviously  the  surface  po¬ 
tential  relative  to  bulk  solution  must  include  the  quadrupole  term.  Figure  5  shows 
the  atomic  charge  density  and  the  electric  field  along  the  z  axis.  Note  that  the  charge 
density  appears  flat  for  Izj  <  0.8  nm.  The  contribution  from  the  ionic  charge  for  z  > 
0.8  nm  is  not  evident  because  the  charge  on  the  water  molecules  dominates.  The 
electric  field  was  obtained  by  integration  of  the  charge  density  from  -oo  to  position 
z.  The  field  is  flat  with  small  variations  around  zero  in  the  region  \z\  <  0.8  nm.  Near 
the  metal  the  electric  field  undergoes  a  series  of  rapid  oscillations  due  to  the  water 
packing  structure  at  the  interface.  The  small  overall  rise  in  field  from  bulk  to  surface 
is  due  to  the  excess  Na^  charge  in  the  screening  layer.  The  oscillation  near  z  =  -1.615 
nm  is  due  to  layering  at  die  restraining  wall.  This  is  an  unwanted  artifact  of  the 
immersed  electrode  model.  It  does  not  occur  if  the  water  density  is  decreased  to  the 
point  that  a  vapor-liquid  interface  opens  up  at  the  restraining  wall.  Space  prevents 
discussion  of  emersed  electrodes[66,  67]  which  have  a  vapor-liquid  interface[68-70] 
and  a  liquid-solid  interface[71]. 

2M  NaCl  Solution 

Figure  6  shows  the  results  for  a  2M  NaCl  solution  at  30°C.  There  are  62  Na  ions, 
58  Cl  ions  and  1516  water  molecules,  the  image  charge  on  the  electrode  is  -4|e|.  The 
cation  and  anion  concentrations  are  approximately  the  same  for  |z|  <  1.0  nm.  The 
bulk  region  appears  larger  than  in  IM  NaCl  solution  consistent  with  narrower 
screening  zone.  There  are  many  similarities  between  Figures  3  and  6.  The  detailed 
proton  and  water  profiles  for  z  >  1 .4  nm  look  the  same.  The  bottom  panel  in  Figure 
6  shows  the  total  potential  calculated  using  the  atom  method,  and  the  components 
of  the  potential  calculated  by  the  molecule  method.  The  contact  potential  is  about 
-2V,  just  the  same  as  in  the  IM  case.  The  difference  between  IM  and  2M  come 
mainly  from  the  ion  distributions  and  their  direct  interaction  with  the  waters.  Thus 
the  monopole  and  dipole  potentials  are  different,  but  since  as  already  seen  for  IM, 
they  cancel  each  other,  the  main  contribution  to  the  total  potential  in  the  bulk  region 
comes  from  the  quadrupole.  The  water  distributions  for  1 M  and  2M  are  also  similar 
and  so  are  the  potentials  which  come  from  water.  Again  if  the  system  were  truly 
neutral  the  monopole  curve  would  be  flat  all  the  way  to  the  edge  of  the  diffuse  layer. 
It  is  flatter  than  the  IM  case  and  oscillates  about  zero  before  diving  down  to  large 
negative  values  for  z  >  1.2  nm.  A  key  observation  is  the  dependence  of  m+d  on 
position.  In  all  the  simulations  the  m+d  potentials  are  similar  even  though  m  and  d 
are  not.  The  transition  to  monotonic  decrease  starts  near  z  =  1.2  nm  is  another  in- 
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Figure  6.  Top.  Probability  density  distribution  profiles  for  2M  NaCl  solution  at  30°C 
and  -4e  electrode  charge.  Bottom.  Total  electric  potential  and  component  multipole 
potentials.  Total  potential  by  atom  method. 


dicator  of  where  the  diffuse  layer  begins  in  this  simulation.  Including  the  quadrupole 
shifts  the  'core'  region  down  by  -3V  and  brings  the  molecular  calculation  of  potential 
into  correspondence  with  the  atom  potential. 

3M  NaCl  Solution 

Figure  7  displays  the  results  for  3M  NaCl.  There  are  94  Na  ions,  86  Cl  ions  and  1463 
water  molecules  at  100°C,  and  the  charge  on  the  electrode  is  -Sjel  or  -0.092  Cm’^. 
Calculations  at  high  temperature  were  originally  selected  to  improve  the  statistics. 
Subsequently  the  temperature  dependence  in  the  range  30  to  100®C  was  found  to  be 
weak.  The  Na"^  screening  charge  is  concentrated  in  the  peak  at  1.45  nm  less  than 
one  water  diameter  removed  from  where  the  9-3  wall  potential  passes  through  zero 
(2^=  1.615  nm).  The  peak  in  the  cation  distribution  results  from  the  strong  primary 
hydration  shell  of  the  cation.  The  solvent  layer  at  the  electrode  also  effects  the  po¬ 
sition  and  shape  of  the  cation  distribution  near  the  metal.  Note  that  for  Iz]  <  1.0  nm 
the  cation  distribution  is  approximately  flat  at  30  ions  nm'*.  There  is  also  a  small 
peak  on  the  left  hand  side  at  ca.  z  =  -1.4  run,  that  is  not  associated  with  screening 
but  is  likely  due  to  layering  of  the  water  molecules  at  the  restraining  wall.  The 
chloride  probability  distribution  has  no  major  structural  features,  certainly  none  like 
the  Na"^  screening  peak  at  z  =  1.45  nm.  Starting  from  the  metal  on  the  right  side  of 
Figure  1,  the  chloride  ion  distribution  rises  to  a  plateau  for  |z|  <  1.00  nm.  The 
chloride  and  sodium  ion  probabilities  are  sufficiently  similar  across  the  plateau  region 
for  us  to  call  this  the  bulk  zone.  This  3M  NaCl  simulation  has  the  best  statistics,  as 
can  be  seen  by  the  degree  of  local  charge  neutrality  (ion  densities  are  the  same),  and 
veiy  nearly  equal  integrated  densities  shown  in  Figure  2. 

The  metal  chaise  is  twice  that  of  IM  and  2M  creating  a  stronger  surface  electric 
field,  which  results  in  more  oriented  waters  in  the  first  layer.  The  height  of  the  proton 
peaks  either  side  of  the  main  water  peak  suggest  that  some  H  bonds  to  the  bulk  region 
are  broken  and  that  OH  bonds  point  directly  toward  the  metal.  This  electric  field 
effect  is  distinct  from  localization  of  water  on  Pt(lOO)  and  Pt(lll)  surfaces  in  the 
simulations  of  Heinzinger  and  Spohr[27],  and  Berkowitz[28,  29],  In  these  papers 
water  is  localized  on  top  sites  of  the  Pt  surface  due  to  directed  features  in  the 
chemisorptive  potential. 

The  bottom  panel  of  Figure  7  shows  the  potential  calculated  using  the  atom 
method,  and  some  of  the  components  of  the  potential  calculated  by  the  molecule 
method.  The  contact  potential  is  larger  due  to  higher  electrode  charge.  Once  again 
the  importance  of  quadrupole  terms  is  apparant.  The  m  curve  due  to  the  ionic  charge 
hovers  around  zero  for  z  <  1 .4  nm.  The  monopole  potential  is  quite  flat  outside  the 
screening  layer.  The  transition  to  monotonic  decrease  starting  near  z  =  1.4  nm  is  an 
indication  of  where  the  diffuse  layer  begins.  Again  m+d  shows  that  the  dipole  po¬ 
tential  completely  compensates  the  monopole  outside  the  screening  zone.  Including 
the  quadrupole  shifts  the  core  region  by  -3V  and  brings  the  m+d  potential  into  closer 
correspondence  with  the  total  potential  calculated  by  the  atom  method. 
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Figure  7.  Top.  Probability  density  distribution  profiles.  3M  NaCl  solution  at  100  C 
and  -8e  electrode  charge.  Bottom.  Total  electric  potential  and  component  multipole 
potentials.  Total  potential  by  atom  method. 


Conclusions 

In  this  paper  we  briefly  reviewed  our  view  the  status  of  MD  simulations  of 
electrochemical  interfaces  and  proposed  a  crude  road  map  of  what  future  calculations 
could  contribute  to  understanding  technology.  As  an  example  of  current  simulation 
capability  we  discussed  the  stmcture  of  electric  double  layers  for  systems  with  1500 
water  molecules  and  salt  concentrations  from  IM  to  3M  NaCl.  Water  stmcture  near 
charged  metal  walls  is  not  an  artifact  of  small  size.  In  IM  NaCl  the  double  layer  is 
about  1  nm  thick  (about  three  layers  of  water)  while  in  3M  solution  the  screening 
layer  was  narrower  than  a  water  molecule.  Water  layers  at  the  surface  significantly 
affected  the  distribution  of  ions  near  the  metal  creating  features  in  the  probability 
distributions  that  are  not  describable  in  the  Gouy-Chapman-Stera  model. 
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